%% set valid mice:
clear validmice
%validmice{1} = '515';validmice{2} = '618';validmice{3} = '953';validmice{4} = '81';validmice{5} = '86';validmice{6} = '146';validmice{7} = '291';validmice{8} = '178';validmice{9} = '149';  %cre negative bmx connexin double knockout
%validmice{1} = '612';validmice{2} = '613';validmice{3} = '619';validmice{4} = '955';validmice{5} = '85';validmice{6} = '147';validmice{7} = '148';validmice{8} = '293';  %cre positive bmx connexin double knockout
%validmice{1} = '515';validmice{2} = '618';validmice{3} = '953';validmice{4} = '81';validmice{5} = '86';validmice{6} = '146';validmice{7} = '291';validmice{8} = '178';validmice{9} = '149';validmice{10} = '612';validmice{11} = '613';validmice{12} = '619';validmice{13} = '955';validmice{14} = '85';validmice{15} = '147';validmice{16} = '148';validmice{17} = '293'; %all aEC connexin mice


%% set analysis type, parameters, thresholds
numbins = 8;
minpoints = 20;     %minimum number of data points for each bin
minmice = 3;        %minimum number of mice to be present in each bin
prestimGcampThreshold = 6;
NAtype = 'distance';        %'distance' is distance from a "positive" gcamp pixel while 'level' is the total gcamp level
takelog = true;     % set this argument to true if the total gcamp measurement is NOT already the log of the total gcamp
vasomotionType = 'maxTotalArteryIOS';  %'maxMaxArteryIOS','maxTotalArteryIOS','maxSurroundIOS','integratedTotalArteryIOS', 'time2peak','maxIOS','totalIOS','integratedMaxVesselIOS'
alpha = .05;
correctForBaseline = false;            %set to true to use baseline value as a fixed effect in LME modeling
modelgenotype = false;
vesselWidthRange = [15,80];     %in microns
dilationcutoff = 100;        %assume you shouldn't be getting > this level of dilation ever
differentFromControlPthresh = .5;
dotrajectoryfiltering = false;
plotMeanTrajectories = false;        %if true, will bootstrap trajectories and plot the result with 95% confidence interval
numBootSamples = 1000;      %number of bootstrap samples to use; only relevant if you're plotting the mean trajectories, but if this number is high then function will be slow
plotcolor = 'k';
frameTime = .14;        %in seconds
prestimtime = 1.5;
baselineSubtract = true;        %this optional measurement will take the baseline value and subtract it from each trial. Mean baseline should already be zero, but since we're normalizing to the GLOBAL baseline, some trials have pre-dilated or pre-constricted vessels
if baselineSubtract   %need to compute a baseline adjustment for the TOTAL dilation as well; that'll depend on the integration window
    intwindow = (11:40)*frameTime;
end

%structure of random effects term
clear lmetype
c = 0;
if correctForBaseline
    c = c+1;
    lmetype{c} = 'width_baseline';
end
if modelgenotype
    c = c + 1;
    lmetype{c} ='genotype';
end

c = c+1;
lmetype{c} = '(1|mouseID:uniqueMovies)';
temp = [];
for n = 1:length(lmetype)
    temp = [temp,' + ',lmetype{n}];
end
lmetype = temp;
%% select binning
edges_distance = [ -1,1,400,800,1200,1600,2000,2500,inf];       %just set up to match 2p

edges_gcamplevel = [0,linspace(2,13,numbins)];

y = nan(numbins,1);
ys = nan(numbins,2);
if modelgenotype
    genotypeeffect = [];
    genotypeeffect_ctrl = [];
end

pooledData = inputtable;        %select some data table to look at; allows manipulation of the table without doing anything to the original

if plotMeanTrajectories || dotrajectoryfiltering || baselineSubtract
    switch vasomotionType
        case 'maxTotalArteryIOS'
            pooledtrajectories = vesselTotalIntensity;
        case 'maxMaxArteryIOS'
            pooledtrajectories = maxInten;
        case 'integratedTotalArteryIOS'
            pooledtrajectories = vesselTotalIntensity;
        case 'maxSurroundIOS'
            pooledtrajectories = surroundInten;
        case 'totalIOS'
            pooledtrajectories = lowResIOSsmooth;
        case 'maxIOS'
            pooledtrajectories = lowResIOSsmooth;
        case 'integratedMaxVesselIOS'
            pooledtrajectories = maxInten;
        case 'time2peak'
            pooledtrajectories = vesselTotalIntensity;
    end

else
    pooledtrajectories = nan(size(pooledData,1),1);
end
if exist('time2peak')
    t2p = time2peak;
elseif ~isempty(strmatch('time2peak',pooledData.Properties.VariableNames))
    t2p = pooledData.time2peak;
else
    t2p = nan(size(pooledtrajectories,1),1);
    t2p(t2p<1) = nan;       %exclude trajectories that have too-short time 2 peak
    t2p(pooledData.evokedMaxTotalVessel<0)=nan;
end

%delete unclassified vessels
t2p = t2p(~isnan(pooledData.vesselID));
pooledtrajectories = pooledtrajectories(~isnan(pooledData.vesselID),:);
pooledData = pooledData(~isnan(pooledData.vesselID),:);
%delete untracked trajectories:
t2p = t2p(~isnan(pooledData.evokedMaxVessel));
pooledtrajectories = pooledtrajectories(~isnan(pooledData.evokedMaxVessel),:);
pooledData = pooledData(~isnan(pooledData.evokedMaxVessel),:);

%take log of total activity (if applicable)
if takelog  
    pooledData.totalgcamp = pooledData.totalgcamp*16;
    pooledData.totalgcamp(pooledData.totalgcamp<5) = 5;     %5 is the activity floor used in NA_Dilation_multilevel_v4_2 so can use this for consistency
    pooledData.totalgcamp = log(pooledData.totalgcamp);
end

%% walk through the table and determine which vessels are differentially responding to local gcamp as compared to the control stimulus
pooledData.differentFromControl = nan(size(pooledData,1),1);
[vesselSegs,~,ic] = unique(pooledData(:,[1,17]),'rows','stable');
pooledData.vesselSegmentID = categorical(ic);       %unique segment ID for each mouse-vessel segment pair
maxControlGcamp = 6;
minEvokedGcamp = 8.5;
minnumpoints = 5;
pvals = nan(size(vesselSegs,1),1);
for n = 1:size(vesselSegs,1)
    cur = pooledData(ic==n,:);
    controltrials = cur(cur.totalgcamp<maxControlGcamp & cur.stimSize == 0,:);
    testtrials = cur(cur.totalgcamp>minEvokedGcamp & cur.stimSize > 0,:);
    %if you want to prepare data for a kruskal-wallis test:
    if sum(~isnan(testtrials.evokedIntegratedTotalVessel)) < minnumpoints || sum(~isnan(controltrials.evokedIntegratedTotalVessel)) < minnumpoints
        continue;
    end
    %if you want to use a rank sum test
    pvals(n) = ranksum(testtrials.evokedIntegratedTotalVessel,controltrials.evokedIntegratedTotalVessel);
    pooledData.differentFromControl(ic==n) = pvals(n);
end

%% run baseline subtraction if wanted:
if baselineSubtract
    bl = mean(pooledtrajectories(:,2:11),2,'omitnan');
    switch vasomotionType
        case 'maxTotalArteryIOS'
            pooledData.evokedMaxTotalVessel = pooledData.evokedMaxTotalVessel-bl;
        case 'maxMaxArteryIOS'
            pooledData.evokedMaxVessel = pooledData.evokedMaxVessel-bl;
        case 'integratedTotalArteryIOS'
            pooledData.evokedIntegratedTotalVessel = pooledData.evokedIntegratedTotalVessel - bl*range(intwindow);
        case 'integratedMaxVesselIOS'
            pooledData.evokedMaxVessel = pooledData.evokedMaxVessel- bl*range(intwindow);
    end
    pooledtrajectories = pooledtrajectories-repmat(bl,1,size(pooledtrajectories,2));
end

%% an option to run additional filtering step on trajectories as is done with the 2p data
if dotrajectoryfiltering
    spectrogramproperties.fs = 1/.14;
    spectrogramproperties.windowsize = 5;
    spectrogramproperties.overlap = 0;
    spectrogramproperties.storedbins = 1:5;
    spectrogramproperties.numfreqpoints = 29;
    stimpowerthresh = 1;

    tempM = mean(pooledtrajectories(:,2:10),2,'omitnan');
    tempS = std(pooledtrajectories(:,2:10),0,2,'omitnan');
    curdddz = pooledtrajectories-repmat(tempM,1,size(pooledtrajectories,2));
    curdddz = curdddz./repmat(tempS,1,size(pooledtrajectories,2));
    [responsive_idx,stimpower,baselinepower,noisybaselinetraj] = trajectoryfilter_v1_1(pooledtrajectories,2:10,11:40,[],[],curdddz,spectrogramproperties);
    responsive_idx = logical(responsive_idx(:,1));      %this value is currently unused
    nondilation = stimpower<stimpowerthresh;
    pooledtrajectories(nondilation,:) = 0;
    pooledData.evokedMaxTotalVessel(~responsive_idx) = nan;
    pooledData.evokedMaxTotalVessel(nondilation) = 0;
    pooledData.maxios(nondilation) = 0;pooledData.totalios(nondilation) = 0;
    pooledData.maxios(~responsive_idx) = nan;pooledData.totalios(~responsive_idx) = nan;
    pooledtrajectories(~responsive_idx,:) = nan;
    pooledtrajectories(nondilation,:) = 0;
    nantrajs = isnan(pooledtrajectories(:,2));
else
    nantrajs = false(size(pooledData,1),1);
end


%% set valid mice:
VM = false(size(pooledData,1),1);
if iscell(validmice) 
    for n = 1:length(validmice)
        try
            VM = VM | contains(pooledData.mouseID,validmice{n});
        catch
            VM = VM | pooledData.mouseID == validmice{n};
        end
    end
else
    VM = validmice;
end
%% some basic housekeeping
% explicitly set categorical variables:
pooledData.mouseID = categorical(pooledData.mouseID);


%get some classifications:
[uniqueMovies,~,uniqueMoviesIdx] = unique(pooledData(:,[1,2,3]),'rows','stable');     %unique movie selection
[uniqueStims,~,uniqueStimsIdx] = unique(pooledData(:,[4,5,6]),'rows','stable');         %unique stimulus selection
[uniqueVessels,~,uniqueVesselsIdx] = unique(pooledData(:,[1,7]),'rows','stable');       %unique vessels selection
pooledData.uniqueMovies = uniqueMoviesIdx;uniqueMovieColumn = size(pooledData,2);                   %include an index identifier for these since you'll need them for the subsequent bootstrapping
pooledData.uniqueStims = uniqueStimsIdx;uniqueStimColumn = size(pooledData,2);
pooledData.uniqueVessels = uniqueVesselsIdx;uniqueVesselColumn = size(pooledData,2);
i_width = pooledData.width_baseline > vesselWidthRange(1) & pooledData.width_baseline < vesselWidthRange(2);
%adjust baseline values so their mean is zero if you want to use them in
%the LME modeling as a fixed effect:
pooledData.vesselTtotalIntensity_baseline = pooledData.vesselTtotalIntensity_baseline - mean(pooledData.vesselTtotalIntensity_baseline,'omitnan');
pooledData.surroundInten_baseline = pooledData.surroundInten_baseline - mean(pooledData.surroundInten_baseline,'omitnan');
pooledData.maxInten_baseline = pooledData.maxInten_baseline - mean(pooledData.maxInten_baseline,'omitnan');
pooledData.width_baseline = pooledData.width_baseline - mean(pooledData.width_baseline,'omitnan');
if ~ismember('genotype',pooledData.Properties.VariableNames)
    modelgenotype = false;
end
%% run through each bin and determine vasomotion parameters:
i_pval = pooledData.differentFromControl < differentFromControlPthresh;

if ismember('stimSize',pooledData.Properties.VariableNames)
    i0 = round(pooledData.stimSize,2) == .23;       %.17 OR .35 OR .52 OR .31(connexin) OR .23 (fine grain receptive field) OR 2.09 for full screen
else
    i0 = pooledData.stimClass == 1;     %for connexin animals, ryr2 animals
end

% I've run likelihood ratio tests comparing including or excluding
% different random effects terms, fixed terms, and baseline width/value is
% definitely important to include, but a unique random term for the
% baseline does not add much

for n = 2:numbins+1
    switch NAtype
        case 'distance'
            xbins = edges_distance;
            i1 = (pooledData.distance2gcamp > edges_distance(n-1)) & (pooledData.distance2gcamp < edges_distance(n));
        case 'level'
            xbins = edges_gcamplevel;
            i1 = (pooledData.totalgcamp > edges_gcamplevel(n-1)) & (pooledData.totalgcamp < edges_gcamplevel(n));
    end
    i2 = pooledData.prestimgcamp < prestimGcampThreshold;
    i3 = pooledData.evokedMaxVessel < dilationcutoff;
    idx = i1 & i2 & i3 & VM & i0 & i_pval & i_width & ~nantrajs;
    curtable = pooledData(idx,:);
    curtraj = pooledtrajectories(idx,:);
    curt2p = t2p(idx);
    it2p = ~isnan(curt2p);
    if size(curtable,1) < minpoints
        continue;
    end
    temp = unique(curtable.mouseID);
    if size(temp,1) < minmice
        continue;
    end
    if plotMeanTrajectories
        % 2 level heirarchical bootstrap:
        if false
            [~,~,ic] = unique(curtable(:,[17,26]),'rows');
            curtable.uniquevessel_stim = categorical(ic);
            IDcolumns = [1,28];     %this should be mouse and roi-stim clusters
            idx = nestedBootstrapNVCsample_v2_0(curtable,numBootSamples,IDcolumns);
            temp = zeros(numBootSamples,size(curtraj,2));
            for n2 = 1:numBootSamples
                for n3 = 1:size(idx,2)
                    temp3 = idx(:,n3,n2);temp3 = temp3(~isnan(temp3));
                    temp(n2,:) = temp(n2,:) + mean(curtraj(temp3,:),'omitnan')/size(idx,2);
                end
            end
            figure
            t  = (1:size(temp,2))*frameTime;t = t-prestimtime;
            cury = median(temp,1);
            curer = [cury-prctile(temp,2.5,1);prctile(temp,97.5,1)-cury];
            shadedErrorBar2(t,cury',curer',plotcolor);
            title([num2str(xbins(n-1)),' < range < ',num2str(xbins(n))])
            beautifyaxis(gcf);
            ylim([-6,30])
            pause(.1)
        elseif true  %3 level heirarchical bootstrap
            
            [~,~,ic] = unique(curtable(:,[uniqueStimColumn,uniqueVesselColumn]),'rows');
            curtable.uniquevessel_stim = categorical(ic);
            [~,~,ic] = unique(curtable(:,[17,uniqueVesselColumn+1]),'rows');
            curtable.uniquevessel_stim_seg = categorical(ic);

            IDcolumns = [1,uniqueVesselColumn+1,uniqueVesselColumn+2];     %this should be mouse and roi-stim clusters
            idx = nestedBootstrapNVCsample_3level_v2_0(curtable,numBootSamples,IDcolumns);
            numframes = size(curtraj,2);
            level1traj = zeros(numBootSamples,size(curtraj,2));
            for n2 = 1:numBootSamples
                level2 = ~isnan(idx(1,1,:,n2));
                level2 = sum(level2(:));
                level2traj = zeros(1,numframes);
                for n3 = 1:level2
                    level3 = ~isnan(idx(1,:,n3,n2));
                    level3 = sum(level3(:));
                    level3traj = zeros(1,numframes);
                    for n4 = 1:level3
                        temp3 = idx(:,n4,n3,n2);temp3 = temp3(~isnan(temp3));
                        level3traj = level3traj + mean(curtraj(temp3,:),'omitnan')/level3;
                    end
                    level2traj = level2traj + level3traj/level2;
                end
                level1traj(n2,:) = level2traj;
            end

            figure
            t  = (1:size(level1traj,2))*frameTime;t = t-prestimtime;
            cury = median(level1traj,1);
            curer = [cury-prctile(level1traj,2.5,1);prctile(level1traj,97.5,1)-cury];
            shadedErrorBar2(t,cury',curer',plotcolor);
            title([num2str(xbins(n-1)),' < range < ',num2str(xbins(n))])
            beautifyaxis(gcf);
            ylim([-6,30])
            pause(.1)
        end
    end
    %so as not to muck with the indexing for the bootstrapping, can append
    %a time to peak column to the table which will show the time it takes
    %to reach the peak dilation
    curtable.time2peak = curt2p;
    switch vasomotionType
        case 'maxMaxArteryIOS'
            lme = fitlme(curtable,['evokedMaxVessel ~ 1',lmetype],'FitMethod','REML');
        case 'maxTotalArteryIOS'
            lme = fitlme(curtable,['evokedMaxTotalVessel ~ 1',lmetype],'FitMethod','REML');
        case 'maxSurroundIOS'
            lme = fitlme(curtable,['evokedMaxSurround ~ 1',lmetype],'FitMethod','REML');
        case 'integratedTotalArteryIOS'
            lme = fitlme(curtable,['evokedIntegratedTotalVessel ~ 1',lmetype],'FitMethod','REML');
        case 'integratedMaxVesselIOS'
            lme = fitlme(curtable,['evokedIntegratedMaxVessel ~ 1',lmetype],'FitMethod','REML');
        case 'maxIOS'
            lme = fitlme(curtable,['maxios ~ 1',lmetype],'FitMethod','REML');
        case 'totalIOS'
            lme = fitlme(curtable,['totalios ~ 1',lmetype],'FitMethod','REML');
        case 'time2peak'
            lme = fitlme(curtable,['time2peak ~ 1',lmetype],'FitMethod','REML');
    end
    [~,~,STATS] = fixedEffects(lme);
    y(n-1) = STATS{1,2};
    feci = coefCI(lme,'Alpha',alpha);
    ys(n-1,1) = y(n-1) - feci(1,1);
    ys(n-1,2) = feci(1,2) - y(n-1);
    if modelgenotype
        if size(STATS,1) < 3
            genotypeeffect = [genotypeeffect;[STATS{2,2},STATS{2,6},abs(STATS{2,2}-STATS{2,7})]];
        else
            genotypeeffect = [genotypeeffect;[STATS{3,2},STATS{3,6},abs(STATS{3,2}-STATS{3,7})]];
        end
    end
end

%% run same analysis but only for control stimulus movies:

if ismember('stimSize',pooledData.Properties.VariableNames)
    i1 = pooledData.stimSize == 0;
else
    disp('missing stim size variable in table')
end

i2 = pooledData.prestimgcamp < prestimGcampThreshold;
i3 = pooledData.evokedMaxVessel < dilationcutoff;
i4 = pooledData.totalgcamp < prestimGcampThreshold;
idx = i1 & i2 & i3 & VM & i4 & i_width;        %with a stimulus window gcamp cutoff PLUS excluding vessels not  different from control

curtable = pooledData(idx,:);
curtraj = pooledtrajectories(idx,:);
curt2p = t2p(idx);
it2p = ~isnan(curt2p);
if plotMeanTrajectories
    % 2 level heirarchical bootstrap:
    if false
        [~,~,ic] = unique(curtable(:,[17,26]),'rows');
        curtable.uniquevessel_stim = categorical(ic);
        IDcolumns = [1,28];     %this should be mouse and roi-stim clusters
        idx = nestedBootstrapNVCsample_v2_0(curtable,numBootSamples,IDcolumns);
        temp = zeros(numBootSamples,size(curtraj,2));
        for n2 = 1:numBootSamples
            for n3 = 1:size(idx,2)
                temp3 = idx(:,n3,n2);temp3 = temp3(~isnan(temp3));
                temp(n2,:) = temp(n2,:) + mean(curtraj(temp3,:),'omitnan')/size(idx,2);
            end
        end
        figure
        t  = (1:size(temp,2))*frameTime;
        cury = median(temp,1);
        curer = [cury-prctile(temp,2.5,1);prctile(temp,97.5,1)-cury];
        shadedErrorBar2(t,cury',curer','g');
        title([num2str(xbins(n-1)),' < range < ',num2str(xbins(n))])
        beautifyaxis(gcf);
        ylim([-6,30])
        pause(.1)
    elseif true  %3 level heirarchical bootstrap
        [~,~,ic] = unique(curtable(:,[uniqueStimColumn,uniqueVesselColumn]),'rows');
        curtable.uniquevessel_stim = categorical(ic);
        [~,~,ic] = unique(curtable(:,[17,uniqueVesselColumn+1]),'rows');
        curtable.uniquevessel_stim_seg = categorical(ic);
        IDcolumns = [1,uniqueVesselColumn+1,uniqueVesselColumn+2];     %this should be mouse and roi-stim clusters

        idx = nestedBootstrapNVCsample_3level_v2_0(curtable,numBootSamples,IDcolumns);
        numframes = size(curtraj,2);
        level1traj = zeros(numBootSamples,size(curtraj,2));
        for n2 = 1:numBootSamples
            level2 = ~isnan(idx(1,1,:,n2));
            level2 = sum(level2(:));
            level2traj = zeros(1,numframes);
            for n3 = 1:level2
                level3 = ~isnan(idx(1,:,n3,n2));
                level3 = sum(level3(:));
                level3traj = zeros(1,numframes);
                for n4 = 1:level3
                    temp3 = idx(:,n4,n3,n2);temp3 = temp3(~isnan(temp3));
                    level3traj = level3traj + mean(curtraj(temp3,:),'omitnan')/level3;
                end
                level2traj = level2traj + level3traj/level2;
            end
            level1traj(n2,:) = level2traj;
        end

        figure
        t  = (1:size(level1traj,2))*.1;
        cury = median(level1traj,1,'omitnan');
        curer = [cury-prctile(level1traj,2.5,1);prctile(level1traj,97.5,1)-cury];
        shadedErrorBar2(t,cury',curer',plotcolor);
        title('control')
        beautifyaxis(gcf);
        ylim([-6,30])
        pause(.1)
    end
end
curtable.time2peak = curt2p;
switch vasomotionType
    
    case 'maxMaxArteryIOS'
        lme = fitlme(curtable,['evokedMaxVessel ~ 1',lmetype],'FitMethod','REML');
    case 'maxTotalArteryIOS'
        lme = fitlme(curtable,['evokedMaxTotalVessel ~ 1',lmetype],'FitMethod','REML');
    case 'integratedMaxVesselIOS'
        lme = fitlme(curtable,['evokedIntegratedMaxVessel ~ 1',lmetype],'FitMethod','REML');
    case 'maxSurroundIOS'
        lme = fitlme(curtable,['evokedMaxSurround ~ 1',lmetype],'FitMethod','REML');
    case 'integratedTotalArteryIOS'
        lme = fitlme(curtable,['evokedIntegratedTotalVessel ~ 1',lmetype],'FitMethod','REML');
    case 'maxIOS'
        lme = fitlme(curtable,['maxios ~ 1',lmetype],'FitMethod','REML');
    case 'totalIOS'
        lme = fitlme(curtable,['totalios ~ 1',lmetype],'FitMethod','REML');
    case 'time2peak'
        curtable = curtable(it2p,:);
        lme = fitlme(curtable,['time2peak ~ 1',lmetype],'FitMethod','REML');
end
[~,~,STATS] = fixedEffects(lme);
yctrl = STATS{1,2};
feci = coefCI(lme,'Alpha',alpha);
ysctrl = [yctrl - feci(1,1);feci(1,2) - yctrl];
if modelgenotype
    if size(STATS,1) < 3
        genotypeeffect_ctrl = [genotypeeffect_ctrl;[STATS{2,2},STATS{2,6},abs(STATS{2,2}-STATS{2,7})]];
    else
        genotypeeffect_ctrl = [genotypeeffect_ctrl;[STATS{3,2},STATS{3,6},abs(STATS{3,2}-STATS{3,7})]];
    end
end

%% plot the results
switch NAtype
    case 'distance'
        plotx = edges_distance(1:(length(edges_distance)-1))+(diff(edges_distance)/2);
    case 'level'
        plotx = edges_gcamplevel(1:(length(edges_gcamplevel)-1))+(diff(edges_gcamplevel)/2);
end
figure;
plotx(end) = plotx(length(plotx)-1)-plotx(length(plotx)-2) + plotx(length(plotx)-1);
shadedErrorBar2(plotx',y,ys(:,1),plotcolor);
hold on
switch NAtype
    case 'distance'
        temp = 2*plotx(end)-plotx(length(plotx)-1);
    case 'level'
        temp = 0;
        xlim([-0.5,14])
end
errorbar(temp,yctrl,ysctrl(1),[plotcolor,'o']);
title(vasomotionType)

if modelgenotype
    plotrange = 1:6;
    figure
    errorbar(plotx(plotrange),-genotypeeffect(plotrange,1),genotypeeffect(plotrange,3),'ko-','linewidth',2,'markerfacecolor','auto')
    hold on
    errorbar(temp,-genotypeeffect_ctrl(1,1),genotypeeffect_ctrl(1,3),'ks','linewidth',2,'markerfacecolor','auto')
    plot([plotx(1),plotx(end)],[0,0],'k--')
    hold off
end
